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Abstract 

Background: Resistance to stress is often heterogeneous among individuals within a population, which helps 
protect against intermittent stress (bet hedging). This is also the case for heat shock resistance in the budding yeast 
Socchoromyces cerevisioe. Interestingly, the resistance appears to be continuously distributed (vs. binary, switch-like) and 
correlated with replicative age (vs. random). Older, slower-growing cells are more resistant than younger, faster-growing 
ones. Is there a fitness benefit to age-correlated stress resistance? 

Results: Here this hypothesis is explored using a simple agent-based model, which simulates a population of individual 
cells that grow and replicate. Cells age by accumulating damage, which lowers their growth rate. They synthesize 
trehalose at a metabolic cost, which helps protect against heat shock. Proteins Tsl 1 and Tps3 (trehalose synthase 
complex regulatory subunit TSL1 and TPS3) represent the trehalose synthesis complex and they are expressed 
using constant, age-dependent and stochastic terms. The model was constrained by calibration and comparison 
to data from the literature, including individual-based observations obtained using high-throughput microscopy 
and flow cytometry. A heterogeneity network was developed, which highlights the predominant sources and 
pathways of resistance heterogeneity. To determine the best trehalose synthesis strategy, model strains with different 
Tsl 1ATps3 expression parameters were placed in competition in an environment with intermittent heat shocks. 

Conclusions: For high severities and low frequencies of heat shock, the winning strain used an age-dependent bet 
hedging strategy, which shows that there can be a benefit to age-correlated stress resistance. The study also illustrates 
the utility of combining individual-based observations and modeling to understand mechanisms underlying population 
heterogeneity, and the effect on fitness. 
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Background 

There is increasing appreciation for individuality of mi- 
crobes [1,2]. Even populations grown up from a single cell, 
in a constant environment can exhibit significant pheno- 
typic heterogeneity in gene expression, protein content and 
physiology. Individual heterogeneity can be important to 
population fitness by allowing different functions (e.g. C 
and N fixation in filamentous cyanobacteria) and survival 
in a fluctuating environment. One prominent example is 
bacterial persistence, where a typical population contains a 
small fraction of slow- or non-growing "persister" cells that 
are not killed by antibiotics [3,4], Cells switch between nor- 
mal and persister states in a random, binary (switch-like) 
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manner. The ability to resist stress comes at a cost (typically 
reduced rates of growth or reproduction). For intermittent 
stress, there is an advantage to maintaining heterogeneity 
among individuals in a population in terms of tradeoffs be- 
tween performance and survival (i.e. an insurance mechan- 
ism referred to as bet hedging [3,5]). 

For eukaryotic microbes, the budding yeast Saccharomy- 
ces cerevisiae (S. cervisiae) is a model organism for stud- 
ying individual heterogeneity and aging and longevity 
[1,6]. Various mechanisms, including stochastic variabil- 
ity in regulatory pathways and production/destruction 
of mRNAs, and deterministic asynchronicity in cell cycle or 
replicative age, lead to heterogeneity in protein content and 
stress resistance in clonal populations [1,7-12], For ex- 
ample, copper resistance is heterogeneous and related 
to cell cycle and replicative age [13]. Intrinsic and in- 
duced expression of heat shock proteins and resistance 
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is heterogeneous [11,14]. Expression of Tsll, used in the 
synthesis of trehalose (an alpha-linked disaccharide of glu- 
cose that enhances thermotolerance, reduces aggregation 
of denatured proteins and protects against oxidative dam- 
age), and heat shock resistance are heterogeneous and 
correlated with replicative age [15]. Natural yeast (i.e. not 
S. cervisiae) populations were found to have heteroge- 
neous resistance to copper, lead and sulfur dioxide, and 
this phenotypic heterogeneity is a beneficial and evolvable 
trait [16]. Unlike in bacteria, heat shock resistance in yeast 
is graded, continuous (vs. binary) and correlated with rep- 
licative age (vs. random). 

The benefit of random heterogeneous expression of a 
stress-response factor has been demonstrated experi- 
mentally and computationally for S. cervisiae [8]. Levy 
et al. [15] hypothesized that correlating heat shock re- 
sistance with age provides an added benefit. The idea is 
that older, slower-growing cells are better candidates for 
being stress resistant because they contribute relatively 
less to the population growth. Is there a fitness benefit 
to age-correlated (vs. random) stress resistance? 

This hypothesis is explored here using modeling, an 
approach that has been applied previously to explore 
the role of heterogeneity [5,17]. The general strategy 
is to develop a mathematical model that includes the 
relevant mechanisms, and then perform numerical 
competition experiments to see if the age-correlated 
resistance trait is beneficial. Such competition and/or 
evolutionary optimization simulations have been used pre- 
viously to determine optimal traits/parameters [18-20]. 
There are many potential model formulations and associ- 
ated parameter sets for simulating replication, aging, re- 
sistance, etc. (reviewed below). To ensure some degree of 
realism, we constrain the model by calibration and com- 
parison to relevant observations from the literature. We 
use agent-based modeling (ABM, aka individual-based 
modeling, IBM), rather than the more common population- 
level modeling approach [18-22]. An ABM is appropriate 
in this case, because it can resolve the continuous/graded 
distribution of various individual properties (e.g., protein 
levels, growth rate, resistance), and model outputs can be 
compared directly to the individual-based observations 
that are used to constrain the model [15,23]. The model 
simulates intracellular mechanisms and the cell behavior 
emerges (systems biology). Then, it simulates many such 
cells and the population behavior emerges (systems ecol- 
ogy). This multi-level approach has been referred to as 
systems bioecology [19-21]. 

We describe the model and compare it to data from 
the literature, which shows that it is generally consistent 
with the observed patterns. A heterogeneity network is 
developed, which highlights the predominant sources 
and pathways of resistance heterogeneity. Then we per- 
form competition experiments with strains that have 



different Tsll/Tps3 expression strategies in an environ- 
ment with intermittent heat shocks. For conditions with 
high severity and low frequency of heat shocks, an age- 
dependent bet hedging strategy is most beneficial, which 
supports the hypothesis of a fitness benefit of age- 
correlated stress resistance. 

Methods 

Model overview 

The model is relatively simple and resolves only those 
mechanisms necessary for exploring the hypothesis and 
comparison to the relevant data. Yeast cells take up glu- 
cose (G, g L" 1 ) and convert it to biomass (Figure 1). Three 
forms of biomass are considered, including structural (m x > 
g dry cell" 1 ), damaged (m D) g dry cell" 1 ) and trehalose (m T , 
g dry cell" 1 ). The total biomass (m, g dry cell" 1 ) is the sum 
of these components (m = rn x + rn D + rn T ). Structural bio- 
mass becomes damaged. A fraction of biomass is synthe- 
sized as trehalose. The model tracks the age or number of 
divisions in terms of bud scars (n B ). A population of indi- 
vidual cells is simulated using an agent-based approach. 

Biomass growth 

A number of metabolism models for S. cervisiae have been 
developed ranging from simple Monod-type growth equa- 
tions to more detailed kinetic models that resolve intracel- 
lular mechanisms up to dynamic/kinetic implementations 
genome-scale network reconstructions [24-28] . 

Here, a relatively simple approach is used, with the 
mass balances for the three components: 




Figure 1 Model schematic. Symbols: G = glucose, m x = structural 
mass, m D = damaged mass, m T = trehalose mass and n B = number 
of bud scars. 
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where ^ (day 1 ) is the growth rate,^ is the trehalose syn- 
thesis fraction and k d (day 1 ) is the damage rate. The 
growth rate increases with glucose and decreases with 
damage: 



«7 



g K g + G K n d d + (m D /m) nd 



(2) 



where ^ mg (day 1 ) is the maximum growth rate, K g (g I/ 1 ) 
is the half-saturation constant for glucose, K d is the half- 
saturation constant and n d is the exponent for damage. 
The model dynamically simulates the extracellular glucose 
concentration considering inflow and washout (for con- 
tinuous culture simulations) and uptake by the cells (see 
Additional file 1: SI text for details). 

Replication 

Cell division in S. cervisiae is via the asymmetrical bud- 
ding process, where a larger mother cell gives birth to a 
smaller daughter cell. With subsequent births, the 
mother s size increases and it accumulates bud scars and 
damage (see below). A number of cell cycle and replica- 
tion models for S. cervisiae have been developed [29,30]. 

Here, the model of Vanoni et al. [30] is adopted. 
Briefly, two cell cycle phases are simulated, unbudded 
and budding. Budding starts at a threshold budding size 
(rn b ), which increases from a specified daughter cell 
value (m b)0 ) by a factor {a m b m nB ~ l ) with each generation. 
Division occurs at a threshold replication size (m r ), 
which is proportional to the budding size (m r =f m>r wib)- 
The daughter gets the mass synthesized during the bud- 
ding phase. The mother gains a bud scar and preferentially 
retains the damage (see below). Individual phenotypic het- 
erogeneity is introduced by randomizing the budding size. 
At birth, the daughter s budding size is drawn from a glo- 
bal truncated normal distribution with specified mean and 
coefficient of variation (CV) [31]. A truncated distribution 
is used to prevent unrealistic values (e.g., m b)0 <0). This 
process introduces non-heritable phenotypic heterogen- 
eity, so values are drawn from a global distribution (vs. 
one with the mean based on the mother). 

Aging and damage accumulation 

Aging in S. cervisiae is due to a number of mechanisms, 
including accumulation of extrachromosomal DNA cir- 
cles (ERCs) and oxidative damage (e.g., carbonylation) to 
proteins [32]. At division, this damage is preferentially 
retained by the mother cell, although the ability to do so 
diminishes with replicative age [33,34]. In addition, the 
mother cell has higher reactive oxygen species (ROS) 
and protein damage rates and lower damaged protein 
degradation rates [35,36]. Trehalose protects against 



ROS damage to proteins [37]. Several generic models of 
aging have been presented [18,34]. Specifically for S. cer- 
visiae, Hirsch [38] developed a model where cells accu- 
mulate a senescence factor at a constant rate and 
partitions it asymmetrically at division. The growth rate 
decreases with increasing amount of this senescence fac- 
tor. More mechanistic models that explicitly represent 
ROS, the damage reaction with a protein (citrate syn- 
thase) and repair reaction with a heat shock protein 
(Hsp90) have been presented [39]. 

The present model considers the production of damaged 
mass (m D ) from structural mass (m x ) in a first-order man- 
ner, at a damage rate that increases with age (kd = a d nB bd ). 
At division, the damage mass is preferentially retained by 
the mother, based on a split fraction (s d ). 

Trehalose synthesis and Tsl1/Tps3 expression 

Trehalose serves as storage carbohydrate and stress pro- 
tectant [40-42]. Synthesis is highest during the stationary 
phase and in response to stress (incl. heat), but trehalose 
also accumulates under normal, non-stressed conditions 
[42-44]. Trehalose is synthesized by a trimeric protein 
complex made up of Tpsl and Tps2, and interchange- 
able Tps3 or Tsll [45,46]. Genes involved in trehalose 
synthesis are induced by heat shock [46,47]. In addition, 
the expression is negatively correlated with growth rate 
and has a stochastic component [15,47]. Tsll and Tps3 
promoters share a common regulatory element (stress- 
responsive element, STRE), but their expression can dif- 
fer [9,10,46]. Trehalose (or more generally carbohydrate 
storage) has been included in metabolic models of S. cer- 
visiae [24,25,48]. 

In the present model, the trehalose synthesis complex 
is represented by Tsll and Tps3, which are assumed to 
limit the trehalose synthesis rate (i.e., Tpsl and Tps2 are 
assumed to be present in excess). Tsll and Tps3 are con- 
sidered separately (vs. the complete synthase complex or 
Tpsl), because their expression differs and to allow for 
direct comparison to observations [15,46]. The trehalose 
synthesis fraction (f s ) is a function of the total trehalose 
synthesis enzyme concentration (e s = e Ts ii + e Tps3 ): 



fs fm 



(3) 



where f m>s is the maximum fraction, K s (mmol I/ 1 ) is the 
half-saturation constant and n s is the exponent for tre- 
halose synthesis. 

Tsll and Tps3 are expressed using a set of constant, age- 
dependent and stochastic terms. In the model applications, 
cells are grown in constant, glucose-replete conditions, so 
the effect of growth condition (i.e., stationary phase) on ex- 
pression is not included. Cells are subjected to heat shock, 
but since the observed resistance was not due to an 
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induced heat shock response [15], induction of resistance 
by heat shock is not included. The Tsll enzyme concentra- 
tion {e Ts i h mol L" 1 ) is assumed to adjust rapidly to the dam- 
age according to: 



ersii 



m D /m 
Krsii + m D /m 



fr. 



Tsll 



(4) 



where e CfTs u (mol L" 1 ) is the magnitude of constant ex- 
pression, e a>Ts u (mol L" 1 ) is the magnitude of damage or 
age-dependent expression, K Ts u is a half-saturation con- 
stant and f r ,Tsii is a randomization factor. Note that ex- 
pression is a function of the combined effect of 
constant, age-dependent and stochastic terms, with their 
relative contribution depending on the assigned param- 
eter values. The randomization factor (f n Tsii) is varied by 
drawing from a global truncated normal distribution 
with mean of 1.0 and specified CV, following the same 
approach used for m^o (see above). An equivalent for- 
mulation is used for Tps3. 

Heat shock tolerance and death 

Heat causes denaturation of proteins and there are a 
number of mechanisms that can prevent this. Trehalose 
stabilizes proteins during heat shock [40]. Other factors 
include various heat shock proteins, whose intrinsic (i.e. 
without heat shock) expression is also heterogeneous 
and correlated with heat shock resistance [11,14], 

Mortality by heat shock is simulated using a determin- 
istic approach [8]. Specifically, the applied heat shock se- 
verity (H a , arbitrary units) is compared to the tolerance 
of the cell (H t ), which increases with the trehalose mass 
fraction (m T I m): 



m T /m 



(5) 



where K h is the half-saturation constant for heat shock 
tolerance. K h is the fraction of trehalose required to 
achieve a tolerance of 0.5. When a heat shock is applied, 
all cells with H t < H a die. H a can be adjusted to reflect 
different experimental conditions. 

Additional mechanisms of heterogeneity 

The model includes a number of deterministic sources of 
heterogeneity, like the uneven split of damage among 
mother and daughter at replication. Also, the budding 
mass threshold (m^o) and Tsll and Tps3 expression fac- 
tors (f r> Tsii and f r> T P s3) are varied stochastically, as described 
above. However, in reality there are numerous other 
mechanisms (e.g., stochastic expression of all genes) that 
contribute to heterogeneity in cellular processes [1,7-12]. 
To account for this, the maximum growth rate (fi mg ) and 
damage exponent (n d ) parameters are also randomized 
(following the same approach used for m b>0 ). 



Agent-based modeling of population 

The model simulates individual yeast cells using an 
agent-based approach [18-22]. Each agent stores the cell 
state variables (e.g., rn x , see Figure 1) and those parame- 
ters that are varied at the individual level (e.g., m b)0 ). For 
continuous culture simulations, the model includes sto- 
chastic washout of cells from the reactor (see Additional 
file 1: SI text for details). Differential equations (e.g., 
Eq. 1) are solved using an explicit numerical integration 
method. The model is implemented in the IAM frame- 
work [19,20], and the source code is available from the 
corresponding author. 

For some experiments the model explicitly simulates 
each individual cell. This includes the microcolony ex- 
periments that have ~10 5 colonies of up to -100 cells 
for a total of 10 7 cells (e.g. Additional file 1: Figure SIB). 
However, for liquid cultures with larger populations, in- 
cluding the competition experiments, this is not feasible. 
For example, a 300-mL culture with a cell density of 2.6 
x 10 8 cells ml/ 1 contains 7.9 x 10 10 cells. Using the 
present model, simulating that many cells for 20 days 
would take approximately 15 years of CPU time and 
18 PB of RAM memory. As is common in microbe ABMs, 
for liquid culture, the model simulates super-individuals, 
which are representative of a number of real individuals 
[21]. Minimum/ maximum numbers of agents are speci- 
fied, and when the number of agents drop/rises below/ 
above this, agents are split/combined (see Additional file 1: 
SI text for details). The number of agents, or the com- 
putational resolution, is set sufficiently high so that the 
model produces robust and reproducible results over 
multiple runs with different seed values for the random 
number generator. 

This application is especially challenging from a compu- 
tational perspective, because of the focus on small frac- 
tions of the population. For example, in one experiment 
by Levy et al. [15], 0.1% of the population was sorted out 
using flow cytometry and the growth rate distribution of 
that fraction was computed (Additional file 1: Figure SI I). 
In order for the model to adequately resolve the hetero- 
geneity of such a small fraction of the population, it needs 
to have a very large number of agents. 

Simulations performed 

The model was constrained by calibration and compari- 
son to relevant observations from the literature. Several 
parameters were calibrated within the available literature 
range with the help of an automated optimization rou- 
tine (see Additional file 1: SI text, Table SI, Figures 
S1&2). Model simulations followed the actual experi- 
mental protocols as described in the respective literature 
references, which could be quite involved. For example, 
one experiment by Levy et al. [15] included growing cells in 
liquid suspension, sub-sampling based on Tsll expression, 
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growing again in liquid suspension, randomly sub- 
sampling, growing as microcolonies, and estimating the 
growth rate of each microcolony based on the change in 
colony area (Additional file 1: Figure SI J). The resulting 
calibrated model is then used without any further changes 
to explore the underlying mechanisms and fitness effect of 
heterogeneous, age -correlated heat shock resistance. 

To understand how heterogeneity is produced and 
how it propagates through the population, we developed 
a heterogeneity network (Figure 2F). The nodes in the 
network represent individual state variables (e.g. m D ), 
calculated variables (e.g. and processes (e.g. unequal 
split of damage, node DAM), and the links represent 
causal relationships. For example, DAM causes heterogen- 
eity in m D) which in turn causes heterogeneity in m (via 
mass summation, node m), ^ (via Eq. 2), and e Ts ii and e Tps3 
(via Eq. 4). By turning off the heterogeneity at a node or 
link in the network and examining the resulting reduction 
in heterogeneity at a downstream node, the heterogeneity 
can be mapped onto the network (see Results section). 

To explore the role of the age-correlated resistance trait 
on the fitness of the yeast, numerical (i.e. simulations) 



competition experiments were performed. Cells were 
grown in a glucose-limited chemostat with constant dilu- 
tion rate (D = 0.15 h" 1 ), a set-up similar to the one used in 
a previous experimental study that examined the effect of 
mild heat shock (28 > 36°C) on S. cervisiae growth rate 
and gene expression [47]. The culture was subjected to 
heat shocks at specified heat shock severity (H a ) and fre- 
quency (F h , h" 1 ). A number of model strains with different 
Tsll/Tps3 expression parameters were developed. For 
Tsll, e C)Ts i2 controls constant expression, e a)Ts u controls 
age-dependent expression, and f r ,Tsii,cv controls stochasti- 
city. Tsll and Tps3 are considered separately in the model 
to allow for comparison to data (Figure 2), but their effect 
on trehalose synthesis is identical (Eq. 3) and a strain with 
a high Tsll parameter (e.g., e c>Ts n) behaves the same as a 
strain with an equivalently high Tps3 parameter (i.e., e c> 
Tps3 ). Therefore, the Tsll and Tps3 parameters were varied 
together (e.g., the ratio e c>Tps3 I e C)Tsil is held constant). We 
created 1,000 model strains, with 10 variants for each Tsll 
parameter. Doing a single simulation with all strains was 
not feasible. Therefore, a tournament-style competition was 
used. In each round, 10 strains were randomly selected and 




Figure 2 Model-data comparison and heterogeneity source and pathway analysis. (A) Damage {m D /m) vs. age (n fi ). Data from [33]. (B) Tsll 
expression (e 7s/7 ) distribution of cells. (C) Age (n e ) vs. Tsll expression (e Tsil ). (D) Heat shock survival of various Tsll -sorted fractions. Data from [15]. "a.u." is 
arbitrary units. (E) Distribution of heat shock tolerance for base case and various diagnostic simulations (e.g. "dam" has equal damage partitioning, 
s d = 05). (F) Heterogeneity network. Line weight indicates contribution of node or link to overall heterogeneity in heat shock tolerance (based on 
variance of normalized H tl e.g., panel E). For details of experiments used to generate the data the reader is referred to the source publications. 
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placed in competition. The winner of each match advanced 
to the next round and the process was repeated until one 
strain was left. For each condition (H a , F h ), the optimal 
bet hedging strategy/parameters emerged as the winner. 
The experiment is implicitly constrained by the metabolic 
cost of trehalose synthesis (reduced synthesis of structural 
biomass, Eq. la), which is traded off against the benefit of 
higher heat shock survival. 

Results and discussion 

Calibration and comparison to data 

The model was calibrated to observations from the 
literature with the help of an automatic calibration 
routine. The database is comprised of 15 datasets 
[9,10,15,33,41,43-45,49]. The reader is referred to the 
original publications for experimental protocols and de- 
tails. This database characterizes the relevant features of 
the system, including the distribution of growth rates, 
damage accumulation with age, Tsll and Tps3 expression 
levels, distribution of Tsll expression, age vs. Tsll expres- 
sion, age distribution for all and top 1% of Tsll expressing 
cells, Tsll expression vs. growth rate, growth rate distribu- 
tions for all, top 1% and 0.1% of Tsll expressing cells, tre- 
halose content, trehalose in wild type vs. Tsll knockout, 
survival vs. Tsll expression, and survival vs. growth rate 
for wild type and Tsll knockout strains. Discussion of all 
datasets is provided in the SI and a selected subset are dis- 
cussed here. 

Oxidative protein damage (carbonyl levels) was ob- 
served to increase with age (bud scars), and the model 
reproduces this general pattern (Figure 2A). The obser- 
vations suggest a step-wise increase whereas the model 
exhibits a more linear shape. The reason for the discrep- 
ancy is unclear. The observations are from a single study 
and it would be useful to obtain additional observations 
to confirm the shape. Damage mass increases with age 
due to preferential inheritance by the mother and an in- 
crease of damage rate with age. The expression of Tsll, 
observed with flow cytometry and green fluorescent protein 
(GFP) fused to Tsll, was quite heterogeneous (Figure 2B). 
The modeled Tsll distribution was not as spread out as the 
observations, for example, the data extended into the nega- 
tive range (presumably a measurement error at low Tsll 
levels), whereas the model restricted Tsll to the positive 
values. The heterogeneity in Tsll is a function of the sto- 
chastic component and the amount of Tsll expression that 
is damage correlated (Eq. 4). Applying the bud scar stain 
WGA-TRITC and passing cells through the flow cytometer 
showed that Tsll expression increased with age (Figure 2C). 
The model also shows an increase. However, for an un- 
known reason it over-predicts the age in the 0.03 Tsll bin. 
When Tsll-sorted sub-populations were heat shocked, sur- 
vival correlated positively with Tsll expression (Figure 2D). 
Again, the model reproduces the general increase, but it 



differs in the 0-0.1% Tsll bin. This narrow bin includes the 
fewest cells and is most susceptible to stochastic variability 
(observations and model). Overall, the model does not 
capture all features of the data, but it reproduced the 
main patterns, including an increase of damage with 
age, heterogeneous Tsll expression, and correlation of 
age and survival with Tsll expression. 

Source and pathway of heterogeneity 

We constructed a heterogeneity network (Figure 2F), 
which defines how heterogeneity can be produced and 
propagate through the population. In the present model, 
all heterogeneity originates at replication. The model does 
not consider other sources of heterogeneity, like stochastic 
differentiation at other times (i.e. in between replication 
events) or heterogeneity in the environment. There are 
a number of deterministic and stochastic sources of 
heterogeneity associated with replication (grey nodes in 
Figure 2F). For example, the scarring process produces 
heterogeneity in bud scars (n B ) in a deterministic manner, 
while the expression of Tsll is varied randomly (f r> Tsii)> 
Despite the numerous sources of heterogeneity, replica- 
tion does not completely randomize or "reset" the cells, 
and the model allows for inter-generation memory. For 
example, bud scars and damage - and thus also the 
growth rate - are heritable (Additional file 1: Figures 
S1A&S3). 

Where does the heterogeneity in survival originate? 
Sequentially removing sources and examining the result- 
ing reduction of heterogeneity in heat shock tolerance 
showed that the scarring and unequal division of damage 
processes are the predominant sources (Figure 2E). But 
there are many ways the heterogeneity can go from these 
sources to heat shock tolerance. How does it propagate 
through the network? Systematically eliminating hetero- 
geneity at links and nodes in the network (e.g., use 
population-average e s in Eq. 3) allowed us to map the 
heterogeneity onto the network. This showed that het- 
erogeneity travels along multiple pathways, but predom- 
inantly from scarring to damage to Tps3 expression to 
trehalose and heat shock tolerance, a deterministic path- 
way that leads to age-correlated stress resistance. 

Competition experiments 

The model did not capture all features of the data, but it 
reproduced the major patterns observed in the relevant 
datasets. It was then used as an experimental system to 
explore the hypothesis outlined in the introduction. To 
determine the best Tsll/Tps3 expression strategy we per- 
formed tournament-style competition experiments be- 
tween 1,000 strains with different expression parameters 
(e c ,Tsii> e a>T swfnTsii,cv) in continuous culture with intermit- 
tent heat shocks. Figure 3A shows the results from one 
simulation at intermediate heat shock severity and 
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Figure 3 Numerical competition experiments. (A) Example of 
one simulation at intermediate heat shock severity and frequency 
{H a = 0.7, F h = 0A4 d~ 1 ). Cell density of 10 competitors. Abrupt drops 
in concentrations correspond to heat shocks. (B) Summary of 
tournament-style competitions. Optimal Tsl 1 /Tps3 expression parameters 
for a number of (B1) heat shock severities {H a varies, F h = 0.14 d" 1 ) and (B2) 
frequencies (H a = 0.7, F h varies). Constant: e cjsjl in uM, Age-dependent: 
e a Tsii in uM, Stochastic: f r j- s n 1 10. Symbols are mean +/- one standard 
deviation of ten replicate experiments. 



frequency. One strain clearly outcompeted the others over 
the course of the experiment. 

We performed a number of experiments at various 
heat shock severities (H a ) and frequencies (F h ), ranging 
from no heat stress {H a = F h = 0) to the maximum heat 
stress the yeast can survive (Figure 3B). When no heat 
shocks were applied, the winning strain had no constant 
or age-dependent terms. It did have a stochastic term, 



but this is simply due to the neutrality of the parameter 
when constant and age-dependent terms are zero (i.e., 
then the stochastic parameter does not affect the expres- 
sion, Eq. 4). At lower heat shock severities, the winning 
strategy was to express Tsll/Tps3 in a constant manner 
without heterogeneity (Figure 3B1). It only takes a small 
amount of trehalose to survive these heat shocks and it 
is best to have all cells synthesize this amount. Adding 
heterogenity would result in some cells being killed by 
heat shock, which would not be beneficial, a finding 
consistent with previous studies [8]. Since heterogeneity 
cannot be avoided with age-dependent expression, con- 
stant expression is the better strategy in this case. At 
higher severities, the amount of trehalose required to 
survive the heat shocks becomes larger and a bet hedg- 
ing strategy becomes beneficial. That is, the average 
amount of trehalose is below what is required to survive 
the heat shocks, but it is heterogeneous and some cells 
have sufficient trehalose to survive the heat shocks, and 
this prevents the population from being wiped out. Under 
such conditions, the model predicted that age-dependent 
expression is better. This can be explained, as suggested in 
the Introduction, by the fact that the older cells contribute 
less to the population growth, and eliminating them is less 
detrimental to population growth (Additional file 1: Figure 
Sll). If age-dependent expression is excluded {e a)Ts u = 0), 
the winning strain has constant and stochastic expression 
terms (Additional file 1: Figure S10). The heterogeneity in- 
troduced through the deterministic aging process is sub- 
optimal and it is beneficial to add more via the stochastic 
term. At lower heat shock frequencies, the winning 
strategy was age-dependent bet hedging, whereas at 
higher frequencies constant expression without hetero- 
geneity was better (Figure 3B2). At lower frequencies, 
the growth period in between the heat shocks is rela- 
tively long and reducing the average trehalose produc- 
tion (as is achieved using a bet hedging strategy) is 
beneficial. At higher frequencies, a bet hedging strategy 
is not advantageous, because too many cells are lost 
through the frequent heat shocks. 

These experiments were performed with a model de- 
signed with equations based on our current understand- 
ing of the underlying mechanisms, a parameter set that 
is generally consistent with the literature and main pat- 
terns consistent with observations. However, we cannot 
rule out that there is not another model formulation (i.e. 
different equations) or parameter set that produces an 
equal-quality calibration but a different result or conclu- 
sion about the fitness effect of age-correlated heat shock 
resistance. This is a common problem in model predic- 
tion and has been referred to as "equifinality" [50], and 
it can be addressed to some extent by varying model for- 
mulations [18,34] and/or parameters [29]. The present 
model is computationally very demanding. Nonetheless, 
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we used an automated optimization routine that allows 
for alternate parameter values, and two runs produced 
essentially the same parameter set, which provides some 
additional support for the robustness of our conclusions. 

Conclusion and outlook 

This study explored the ecological role of heterogeneous, 
age-correlated heat shock resistance in S. cervisiae. A sim- 
ple model was constructed based on our current under- 
standing of the underlying mechanisms, and comparison to 
relevant data shows it is consistent with observed patterns. 
Competition experiments with strains that have different 
stress protectant synthesis strategies shows that, for high 
severities and low frequencies of heat shock, an age- 
dependent bet hedging strategy is best. This supports the 
hypothesis that age-correlated resistance is more beneficial 
than random resistance. Although the model is specific to 
heat shock resistance in S. cervisiae, trehalose is produced 
by many different organisms and also protects against other 
forms of stress (e.g., ethanol, [41]), so our results have 
broader relevance. However, there are also cases where 
resistance is negatively correlated with age (i.e., younger 
cells are more resistant), like Sodlp-mediated copper 
resistance [13], so these results cannot be generalized to 
all types of stress. 

The finding that it is advantageous for older cells to in- 
vest in increasing stress tolerance has implication for un- 
derstanding aging and longevity- two very different things, 
with different selective forces acting [6]. Longevity is a 
highly adaptive trait and it is generally considered that 
genes promoting longer lifespan do so by improving som- 
atic resistance in unfavorable conditions [6]. Our results 
provide a clear example of how such a mechanism could 
operate. 

Our model was designed specifically for exploring the 
role of age-correlated heat shock resistance in S. cervisiae. 
For that purpose it was kept as simple as possible, while still 
including the relevant mechanisms. This naturally limits 
the models applicability to other questions, although it 
should be useful for exploring other features related to 
aging, heterogeneity and stress resistance. For example, 
with minimal changes (i.e. Eq. 4), the present model could 
be used to predict expression of other proteins. The model 
can also serve as a stepping stone for further model deve- 
lopment. A lot more is known about the various mech- 
anisms involved in the problem and this knowledge is 
sufficient to support the development of a more de- 
tailed model. It would be interesting to bring in more 
mechanistically-detailed models of gene transcription 
and expression noise [7,8,12], more detailed and/or 
genome-scale metabolism [24,26-28] and cell cycle 
control [29]. Sub-genome scale combined signaling, gene 
expression and metabolism models have been developed 
[51]. It seems that several pieces are in place to support the 



development of such a model, which would require a large 
community effort (as was done for the latest metabolic net- 
work reconstruction, [27]), but it would be worth it. 

This study combined individual-based observations (IBO) 
and modeling (IBM) to understand mechanisms underlying 
population heterogeneity, and the effect on fitness [23]. 
Individual-based observation and experimental tech- 
nologies are advancing rapidly and are generating large 
amount of novel data [15,52]. These data are different 
from traditional population-level observations, which 
were amenable to analysis using traditional population- 
level models, and they require new methods and models. 
Our study illustrates the utility of combining IBM and 
IBO. The IBOs of Levy et al. [15] were used to constrain 
the individual-level processes in the IBM. The IBM, in 
turn, put the IBOs into ecological context. 

This paper presents the use of a multi-scale modeling 
approach to investigate the role of an intracellular mech- 
anism in the ecological fitness of an organism. Covering 
multiple levels of organization is a general problem in 
the biological sciences. Several systems approaches have 
been developed to address this challenge [53,54]. The 
approach used here, "systems bioecology", combines sys- 
tems biology and systems ecology [19-21]. The idea is 
conceptually quite simple. First, the intracellular states 
and mechanisms of microorganisms are explicitly simu- 
lated (systems biology). Then, whole populations of these 
individual microbes are simulated directly using agent- 
based modeling (ABM), including their interaction with 
the environment (systems ecology). This general approach 
may be applicable to other questions involving the role of 
intracellular mechanisms at the ecosystem scale. 

Additional file 



Additional file 1: Additional model details. Additional discussion of 
model results. 
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